In [12]:
using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations, DifferentialEquations
gr();
Dr. Sheehan Olver
s.olver@imperial.ac.uk
Website: https://github.com/dlfivefifty/M3M6LectureNotes
This lecture we do the following:
That is, we introduce recurrences related to ultraspherical polynomials. This allows us to represent general linear differential equations as almost-banded systems.
We have discussed general properties, but now we want to discuss some classical orthogonal polynomials, beginning with Chebyshev (first kind) $T_n(x)$, which is orthogonal w.r.t. $$1\over \sqrt{1-x^2}$$ and ultraspherical $C_n^{(\lambda)}(x)$, which is orthogonal w.r.t. $$(1-x^2)^{\lambda - \half}$$ for $\lambda > 0$. Note that Chebyshev (second kind) satisfies $U_n(x) = C_n^{(1)}(x)$.
For Chebyshev, recall we have the normalization constant (here we use a superscript $T_n(x) = k_n^{\rm T} x^n + O(x^{n-1})$) $$ k_0^{\rm T} = 1, k_n^{\rm T} = 2^{n-1} $$ For Ultraspherical $C_n^{(\lambda)}$, this is $$ k_n^{(\lambda)} = {2^n (\lambda)_n \over n!} = {2^n \lambda (\lambda+1) (\lambda+2) \cdots (\lambda+n-1) \over n!} $$ where $(\lambda)_n$ is the Pochhammer symbol. Note for $U_n(x) = C_n^{(1)}(x)$ this simplifies to $k_n^{\rm U} = k_n^{(1)} = 2^n$.
We have already found the recurrence for Chebyshev: $$ x T_n(x) = {T_{n-1}(x) \over 2} + {T_{n+1}(x) \over 2} $$ We will show that we can use this to find the recurrence for all ultraspherical polynomials. But first we need some special recurrences.
Remark Jacobi, Laguerre, and Hermite all have similar relationships, which will be discussed further in the problem sheet.
It turns out that the derivative of $T_n(x)$ is precisely a multiple of $U_{n-1}(x)$, and similarly the derivative of $C_n^{(\lambda)}$ is a multiple of $C_{n-1}^{(\lambda+1)}$.
Proposition (Chebyshev derivative) $$T_n'(x) = n U_{n-1}(x)$$
Proof We first show that $T_n'(x)$ is othogonal w.r.t. $\sqrt{1-x^2}$ to all polynomials of degree $m < n-1$, denoted $f_m$, using integration by parts: $$ \ip<T_n',f_m>_{\rm U} = \int_{-1}^1 T_n'(x) f_m(x) \sqrt{1-x^2} \dx = -\int_{-1}^1 T_n(x) (f_m'(x)(1-x^2) + xf_m) {1 \over \sqrt{1-x^2}} \dx = - \ip<T_n, f_m'(1-x^2) + x f_m >_{\rm T} = 0 $$ since $f_m'(1-x^2) + f_m $ is degree $m-1 +2 = m+1 < n$.
The constant works out since $$ T_n'(x) = {\D \over \dx} (2^{n-1} x^n) + O(x^{n-2}) = n 2^{n-1} x^{n-1} + O(x^{n-2}) $$ ⬛️
The exact same proof shows the following:
Proposition (Ultraspherical derivative) $${\D \over \dx} C_n^{(\lambda)}(x) = 2 \lambda C_{n-1}^{(\lambda+1)}(x)$$
Like the three-term recurrence and Jacobi operators, it is useful to express this in matrix form. That is, for the derivatives of $T_n(x)$ we get $$ {\D \over \dx} \begin{pmatrix} T_0(x) \\ T_1(x) \\ T_2(x) \\ \vdots \end{pmatrix}= \begin{pmatrix} 0 \cr 1 \cr & 2 \cr && 3 \cr &&&\ddots \end{pmatrix} \begin{pmatrix} U_0(x) \\ U_1(x) \\ U_2(x) \\ \vdots \end{pmatrix} $$ which let's us know that, for $$ f(x) = (T_0(x),T_1(x),\ldots) \begin{pmatrix} f_0\\f_1\\\vdots \end{pmatrix} $$ we have a derivative operator in coefficient space as $$ f'(x) = (U_0(x),U_1(x),\ldots)\begin{pmatrix} 0 & 1 \cr && 2 \cr &&& 3 \cr &&&&\ddots \end{pmatrix} \begin{pmatrix} f_0\\f_1\\\vdots \end{pmatrix} $$
Demonstration Here we see that applying a matrix to a vector of coefficients successfully calculates the derivative:
In [13]:
f = Fun(x -> cos(x^2), Chebyshev()) # f is expanded in Chebyshev coefficients
n = ncoefficients(f) # This is the number of coefficients
D = zeros(n-1,n)
for k=1:n-1
D[k,k+1] = k
end
D
Out[13]:
Here D*f.coefficients gives the vector of coefficients corresponding to the derivative, but now the coefficients are in the $U_n(x)$ basis, that is, Ultraspherical(1):
In [14]:
fp = Fun(Ultraspherical(1), D*f.coefficients)
fp(0.1)
Out[14]:
Indeed, it matches the "true" derivative:
In [15]:
f'(0.1)
Out[15]:
In [16]:
-2*0.1*sin(0.1^2)
Out[16]:
Note that in ApproxFun.jl we can construct these operators rather nicely:
In [17]:
D = Derivative()
(D*f)(0.1)
Out[17]:
Here we see that we can write produce the ∞-dimensional version as follows:
In [18]:
D : Chebyshev() → Ultraspherical(1)
Out[18]:
We can convert between any two polynomial bases using a lower triangular operator, because their span's are equivalent. In the case of Chebyshev and ultraspherical polynomials, they have the added property that they are banded.
Proposition (Chebyshev T-to-U conversion) \begin{align*} T_0(x) &= U_0(x) \\ T_1(x) &= {U_1(x) \over 2} \\ T_n(x) &= {U_n(x) \over 2} - {U_{n-2}(x) \over 2} \end{align*}
Proof
Before we do the proof, note that the fact that there are limited non-zero entries follows immediately: if $m < n-2$ we have $$ \ip<T_n, U_m>_{\rm U} = \ip<T_n, (1-x^2)U_m>_{\rm T} = 0 $$
To actually determine the entries, we use the trigonometric formulae. Recall for $x = (z + z^{-1})/2$, $z = \E^{\I \theta}$, we have \begin{align*} T_n(x) &= \cos n \theta = {z^{-n} + z^n \over 2}\\ U_n(x) &= {\sin (n+1) \theta \over \sin \theta} = {z^{n+1} - z^{-n-1} \over z - z^{-1}} = z^{-n} + z^{2-n} + \cdots + \cdots + z^{n-2} + z^n \end{align*} The result follows immediately.
⬛️
Corollary (Ultrapherical λ-to-(λ+1) conversion) $$ C_n^{(\lambda)}(x) = {\lambda \over n+ \lambda} (C_n^{(\lambda+1)}(x) - C_{n-2}^{(\lambda+1)}(x)) $$
Proof This follows from differentiating the previous result. For example: \begin{align*} {\D\over \dx} T_0(x) &= {\D\over \dx} U_0(x) \\ {\D\over \dx} T_1(x) &= {\D\over \dx} {U_1(x) \over 2} \\ {\D\over \dx} T_n(x) &= {\D\over \dx} \left({U_n(x) \over 2} - {U_{n-2} \over 2} \right) \end{align*} becomes \begin{align*} 0 &= 0\\ U_0(x) &= C_1^{(2)}(x) \\ n U_{n-1}(x) &= C_{n-1}^{(2)}(x) - C_{n-3}^{(2)}(x) \end{align*}
Differentiating this repeatedly completes the proof.
⬛️
Note we can write this in matrix form, for example, we have $$ \underbrace{\begin{pmatrix}1 \cr 0 & \half\cr -\half & 0 & \half \cr &\ddots &\ddots & \ddots\end{pmatrix} }_{S_0^\top} \begin{pmatrix} U_0(x) \\ U_1(x) \\ U_2(x) \\ \vdots \end{pmatrix} = \begin{pmatrix} T_0(x) \\ T_1(x) \\ T_2(x) \\ \vdots \end{pmatrix} $$
therefore, $$ f(x) = (T_0(x),T_1(x),\ldots) \begin{pmatrix} f_0\\f_1\\\vdots \end{pmatrix} = (U_0(x),U_1(x),\ldots) S_0 \begin{pmatrix} f_0\\f_1\\\vdots \end{pmatrix} $$
Again, we can construct this nicely in ApproxFun:
In [19]:
S₀ = I : Chebyshev() → Ultraspherical(1)
Out[19]:
In [20]:
f = Fun(exp, Chebyshev())
g = S₀*f
Out[20]:
In [21]:
g(0.1) - exp(0.1)
Out[21]:
Theorem (three-term recurrence for Chebyshev U) \begin{align*} x U_0(x) &= {U_1(x) \over 2} \\ x U_n(x) &= {U_{n-1}(x) \over 2} + {U_{n+1}(x) \over 2} \end{align*}
Proof Differentiating \begin{align*} x T_0(x) &= T_1(x) \\ x T_n(x) &= {T_{n-1}(x) \over 2} + {T_{n+1}(x) \over 2} \end{align*} we get \begin{align*} T_0(x) &= U_0(x) \\ T_n(x) + n x U_{n-1}(x) &= {(n-1) U_{n-2}(x) \over 2} + {(n+1) U_n(x) \over 2} \end{align*} substituting in the conversion $T_n(x) = (U_n(x) - U_{n-2}(x))/2$ we get \begin{align*} T_0(x) &= U_0(x) \\ n x U_{n-1}(x) &= {(n-1) U_{n-2}(x) \over 2} + {(n+1) U_n(x) \over 2} - (U_n(x) - U_{n-2}(x))/2 = {n U_{n-2}(x) \over 2} + {n U_n(x) \over 2} \end{align*}
⬛️
Differentiating this theorem again and applying the conversion we get the following
Corollary (three-term recurrence for ultrashperical) \begin{align*} x C_0^{(\lambda)}(x) &= {1 \over 2\lambda } C_1^{(\lambda)}(x) \\ x C_n^{(\lambda)}(x) &= {n+2\lambda-1 \over 2(n+\lambda)} C_{n-1}^{(\lambda)}(x) + {n+1 \over 2(n+\lambda)} C_{n+1}^{(\lambda)}(x) \end{align*}
Here's an example of the Jacobi operator (which is the transpose of the multiplciation by $x$ operator):
In [22]:
Multiplication(Fun(), Ultraspherical(2))'
Out[22]:
The preceding results allowed us to represent
as banded operators. We will see that we can combine these, along with
4. Evaluation
to solve ordinary differential equations.
Consider the simplest ODE: \begin{align*} u(0) &= 0 \\ u'(x) - u(x) &= 0 \end{align*} and suppose represent $u(x)$ in its Chebyshev expansion, with to be determined coefficents. In other words, we want to calculate coefficients $u_k$ such that $$ u(x) = \sum_{k=0}^\infty u_k T_k(x) = (T_0(x), T_1(x), \ldots) \begin{pmatrix} u_0 \\ u_1 \\ \vdots \end{pmatrix} $$ In this case we know that $u(x) = \E^x$, but we would still need other means to calculate $u_k$ (They are definitely not as simple as Taylor series coefficients).
We can express the constraints as acting on the coefficients. For example, we have $$ u(0) = (T_0(0), T_1(0), \ldots) \begin{pmatrix} u_0\\u_1\\\vdots \end{pmatrix} = (1,0,-1,0,1,\ldots) \begin{pmatrix} u_0\\u_1\\\vdots \end{pmatrix} $$ We also have $$u'(x) = (U_0(x),U_1(x),\ldots) \begin{pmatrix} 0 & 1 \cr && 2 \cr &&& 3 \cr &&&&\ddots \end{pmatrix}\begin{pmatrix} u_0\\u_1\\\vdots \end{pmatrix} $$ To represent $u'(x) - u(x)$, we need to make sure the bases are compatible. In other words, we want to express $u(x)$ in it's $U_k(x)$ basis using the conversion operator $S_0$: $$u(x) = (U_0(x),U_1(x),\ldots) \begin{pmatrix} 1 &0 & -\half \cr & \half & 0 & -\half \cr &&\ddots & \ddots & \ddots \end{pmatrix}\begin{pmatrix} u_0\\u_1\\\vdots \end{pmatrix} $$
Which gives us, $$ u'(x) - u(x) = (U_0(x),U_1(x),\ldots) \begin{pmatrix} -1 &1 & \half \cr & -\half & 2 & \half \cr && -\half & 3 & \half \cr &&&\ddots & \ddots & \ddots \end{pmatrix} \begin{pmatrix} u_0\\u_1\\\vdots \end{pmatrix} $$
Combing the differential part and the evaluation part, we arrive at an (infinite) system of equations for the coefficients $u_0,u_1,\dots$: $$ \begin{pmatrix} 1 & 0 & -1 & 0 & 1 & \cdots \\ -1 &1 & \half \cr & -\half & 2 & \half \cr && -\half & 3 & \half \cr &&&\ddots & \ddots & \ddots \end{pmatrix} \begin{pmatrix} u_0\\u_1\\\vdots \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \\ 0 \\ \vdots \end{pmatrix} $$
How to solve this system is outside the scope of this course (though a simple approach is to truncate the infinite system to finite systems). We can however do this in ApproxFun:
In [23]:
B = Evaluation(0.0) : Chebyshev()
D = Derivative() : Chebyshev() → Ultraspherical(1)
S₀ = I : Chebyshev() → Ultraspherical(1)
L = [B;
D - S₀]
Out[23]:
We can solve this system as follows:
In [25]:
u = L \ [1; 0]
plot(u)
Out[25]:
It matches the "true" result:
In [26]:
u(0.1) - exp(0.1)
Out[26]:
Note we can incorporate right-hand sides as well, for example, to solve $u'(x) - u(x) = f(x)$, by expanding $f$ in its Chebyshev U series.
This approach extends to second-order constant-coefficient equations by using ultraspherical polynomials. Consider \begin{align*} u(-1) &= 1\\ u(1) &= 0\\ u''(x) + u'(x) + u(x) &= 0 \end{align*} Evaluation works as in the first-order case. To handle second-derivatives, we need $C^{(2)}$ polynomials:
In [27]:
D₀ = Derivative() : Chebyshev() → Ultraspherical(1)
D₁ = Derivative() : Ultraspherical(1) → Ultraspherical(2)
D₁*D₀ # 2 zeros not printed in (1,1) and (1,2) entry
Out[27]:
For the identity operator, we use two conversion operators:
In [28]:
S₀ = I : Chebyshev() → Ultraspherical(1)
S₁ = I : Ultraspherical(1) → Ultraspherical(2)
S₁*S₀
Out[28]:
And for the first derivative, we use a derivative and then a conversion:
In [29]:
S₁*D₀ # or could have been D₁*S₀
Out[29]:
Putting everything together we get:
In [30]:
B₋₁ = Evaluation(-1) : Chebyshev()
B₁ = Evaluation(1) : Chebyshev()
# u(-1)
# u(1)
# u'' + u' +u
L = [B₋₁;
B₁;
D₁*D₀ + S₁*D₀ + S₁*S₀]
Out[30]:
In [31]:
u = L \ [1.0,0.0,0.0]
plot(u)
Out[31]:
In [32]:
x = Fun()
Jᵗ = Multiplication(x) : Chebyshev() → Chebyshev() # transpose of the Jacobi operator
Out[32]:
We set op ther system as follows:
In [33]:
L = [B₋₁; # u(-1)
B₁ ; # u(1)
D₁*D₀ - S₁*S₀*Jᵗ] # u'' - x*u
Out[33]:
In [34]:
u = L \ [1.0;0.0;0.0]
plot(u; legend=false)
Out[34]:
If we introduce a small parameter, that is, solve \begin{align*} u(-1) &= 1\\ u(1) &= 0\\ \epsilon u''(x) - xu(x) &= 0 \end{align*} we can see pretty hard to compute solutions:
In [35]:
ε = 1E-6
L = [B₋₁;
B₁ ;
ε*D₁*D₀ - S₁*S₀*Jᵗ]
u = L \ [1.0;0.0;0.0]
plot(u; legend=false)
Out[35]:
Because of the banded structure, this can be solved fast:
In [37]:
ε = 1E-10
L = [B₋₁;
B₁ ;
ε*D₁*D₀ - S₁*S₀*Jᵗ]
@time u = L \ [1.0;0.0;0.0]
@show ncoefficients(u);
To handle other variable coefficients, first consider a polynomial $p(x)$. If Multiplication by $x$ is represented by multiplying the coefficients by $J^\top$, then multiplication by $p$ is represented by $p(J^\top)$:
In [38]:
M = -I + Jᵗ + (Jᵗ)^2 # represents -1+x+x^2
Out[38]:
In [39]:
ε = 1E-6
L = [B₋₁;
B₁ ;
ε*D₁*D₀ - S₁*S₀*M]
@time u = L \ [1.0;0.0;0.0]
@show ε*u''(0.1) - (-1+0.1+0.1^2)*u(0.1)
plot(u)
Out[39]:
For other smooth functions, we first approximate in a polynomial basis, and without loss of generality we use Chebyshev T basis. For example, consider \begin{align*} u(-1) &= 1\\ u(1) &= 0\\ \epsilon u''(x) - \E^x u(x) &= 0 \end{align*} where $$ \E^x \approx p(x) = \sum_{k=0}^{m-1} p_k T_k(x) $$ Evaluating at a point $x$, recall Clenshaw's algorithm: \begin{align*} \gamma_{n-1} &= 2p_{n-1} \\ \gamma_{n-2} &= 2p_{n-2} + 2x \gamma_{n-1} \\ \gamma_{n-3} &= 2 p_{n-3} + 2x \gamma_{n-2} - \gamma_{n-1} \\ & \vdots \\ \gamma_1 &= p_1 + x \gamma_2 - \half \gamma_3 \\ p(x) = \gamma_0 &= p_0 + x \gamma_1 - \half \gamma_2 \end{align*} If multiplication by $x$ becomes $J^\top$, then multiplication by $p(x)$ becomes $p(J^\top)$, and hence we calculate:\ \begin{align*} \Gamma_{n-1} &= 2p_{n-1}I \\ \Gamma_{n-2} &= 2p_{n-2}I + 2J^\top \Gamma_{n-1} \\ \Gamma_{n-3} &= 2 p_{n-3}I + 2J^\top \Gamma_{n-2} - \Gamma_{n-1} \\ & \vdots \\ \Gamma_1 &= p_1I + J^\top \Gamma_2 - \half \Gamma_3 \\ p(J^\top) = \Gamma_0 &= p_0 + x \Gamma_1 - \half \Gamma_2 \end{align*}
Here is an example:
In [40]:
p = Fun(exp, Chebyshev()) # polynomial approximation to exp(x)
M = Multiplication(p) : Chebyshev() # constructed using Clenshaw:
Out[40]:
In [41]:
ApproxFun.bandwidths(M) # still banded
Out[41]:
In [42]:
ε = 1E-6
L = [B₋₁;
B₁ ;
ε*D₁*D₀ + S₁*S₀*M]
@time u = L \ [1.0;0.0;0.0]
@show ε*u''(0.1) + exp(0.1)*u(0.1)
plot(u)
Out[42]: